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Abstract. For a minimal Hubbard-type system at different interaction strengths U, we 
investigate the density-response for an excitation beyond the linear regime using the generalized 
Kadanoff-Baym ansatz (GKBA) and the second Born (2B) approximation. We find strong 
correlation features in the response spectra and establish the connection to an involved double 
excitation process. By comparing approximate and exact Green's function results, we also 
observe an anomalous {/-dependence of the energy of this double excitation in 2B+GKBA. 
This is in accordance with earlier findings [K. Balzer et ai, EPL 98, 67002 (2012)] on double 
excitations in quantum wells. 



1. Introduction 

The study of strongly correlated quantum systems and materials, e.g., |1J and references therein, 
is a very rapidly evolving field in both experimental and theoretical physics. Especially the out- 
of-equilibrium dynamics is of great current interest in solid-state, atomic and molecular physics, 
in nanoelectronics, quantum transport etc. In all these fields, the availability of intense and 
coherent radiation, combined with ultra-short laser pulses, has triggered many key experiments 
that allow one to investigate matter under extreme nonequilibrium conditions where correlation 
and nonlinear effects are important. As examples, consider the photoionization of multi-electron 
atoms and molecules [2] , the many-body dynamics of particles in optical lattices [3] or quantum 
interference effects in Mott insulators [JJ. 

In the recent decade, the nonequilibrium Green's function (NEGF) approach has developed 
into a powerful numerical tool and is able to address the dynamics of very different quantum 
systems. Examples include the excitation of non-ideal (semiconductor) plasmas [3|6], nuclear 
collisions [7j, the carrier dynamics and carrier-phonon interaction in quantum dots, wells 
and lattice networks contacted to leads [H O [10l [11], and the problem of baryogenesis in 
cosmology [12]. In all these cases, the central quantity is the one-particle nonequilibrium Green's 



function denned on the complex Keldysh contour C, 



g° s ,(t,t') = -i{T c c s , a (t)c\ l>CJ {t')) (1) 
= c (t- <&>(i, + 9 c (t' - t) tf<(t, t>) , 

with non-temporal degrees of freedom s, s' and a. For our purpose, s and s' will refer to spatial 
degrees of freedom, and a =t, i will indicate the spin, assuming diagonality in the spin subspace. 
Further, the operator Tq accounts for contour ordering of the times t and t' , and (. . .) means 
averaging in the grand canonical ensemble. 

In principle, the two-time dependence of the Green's function enables the systematic and 
exact treatment of dynamic correlation effects independently of the strength of the interaction 
and regardless of the presence of weak or strong external fields. Even interaction quenches can 
be studied, e.g. p3]. However, the application to large system sizes and the investigation of the 
system's long time behavior usually demand additional approximations and simplifications to 
the equations of motions of the NEGF |14| . 

The most straightforward approach is to perform a many-body approximation (MBA) on 
the self-energy. Due to diagram expansions, this can be done in a systematic and conserving 
way and leads, e.g., to the Hartree-Fock (HF), second (-order) Born (2B), GW or T-matrix 
approximations. Correlation effects are then partially treated beyond the level of Hartree- 
Fock. Additional simplifications typically concern the two-time structure which renders NEGF 
calculations demanding |15tll6j. Within the generalized Kadanoff Baym ansatz (GKBA) |17tll8|. 
it is possible to study the equations of motion for the Green's function in the single-time limit, 
which drastically reduces the numerical effort. 

In general, the application of MBAs and other simplifications such as the GKBA require 
thorough justification as they can lead to unphysical effects such as self- inter action errors and 
spurious dynamical excitations [19], bistability [11] or artificial steady states [20]. For this 
reason, we, in this contribution, extend previous work on weak perturbations [21], which included 
the discussion of double excitations (see also |22|). to the nonlinear regime and analyze the 
performance of the GKBA for a finite system with the two-body interactions being treated in 
the second Born approximation. 



2. Equations of motion and GKBA 

Introducing the one-particle self-energy Y,° s i(t, t'), the Green's function ([!]) obeys the Kadanoff- 
Baym equation (KBE) [23] . 

" 9-sAt,t') =Sc(t- t')5 ss , + jf dt £^(i,i>?,(M') , (2) 

and its adjoint with t t'. Here, h a i denotes the single-particle (kinetic plus potential) energy, 
and summation over s is implied. 

Using the generalized Kadanoff-Baym ansatz, we transform the set of KBEs into a single 
equation for the reduced density matrix p°' s f(t) = —ig^f^jt) or equivalently p^(t) = 
—ig°'j'(t,t) retaining time causality and important conservation laws such as the conservation 
of total energy, momentum and density. In the course of this, the time-off-diagonal components 
of the NEGF are reconstructed as (we again sum over s) , 

9 a s f(t,t') = - 9 T\t,t>)p^(t') + p^(t)g°f v (t,t>) , (3) 

and inserted into the time-diagonal collision integral (the r.h.s. of the KBEs for t = t' in 
the second Born approximation). For an explicit expression of £"_,(£,£') as functional of the 



NEGF and for the corresponding collision term, see, e.g., the contribution of S. Hermanns et 
al. in the present conference proceedings. Finally, the treatment of the two-time retarded and 
advanced propagators ggg/ Gt (t,tf) and g°' s , v (t,t?) in Hartree-Fock approximation renders the 
ansatz practical: 

(4) 

ss' 

Here, h^ F (t) denotes the single-particle time-dependent Hartree-Fock Hamiltonian. 

We emphasize that the reconstruction of the greater and lesser components of the Green's 
function with real t and t' is sufficient as long as the method of adiabatic switching is applied 
to generate the correlated initial (ground) state by time propagation, for details see [18J. 

3. Nonlinear density-response 

As outlined above, we are interested in the dynamical behavior of a finite quantum system beyond 
the regime of linear response. To this end, we consider a one-dimensional Hubbard model at 
half-filling with hopping T and on-site interaction U and strongly perturb it by changing the 
site energies in time [24] . 

The initial Hamiltonian for times t < reads, 

Ho = ~T fi t<^V + u J2 fl s,t fl s,i , (5) 

<s,s'> cr=t4 s 

where s and s' range from to L — 1 for a chain of length L, and < s, s' > indicates nearest- 
neighbor sites. Further, n SjCT = cl >a c s>rT denotes the density operator, and the energy (time) 
is measured in units of T (the inverse hopping T _1 ). Generally, we study the chain for open 
boundary conditions but this is irrelevant for the special case of two sites. 

The system is excited by an instantaneous change of the site energies from zero to a positive 
finite value e s . If two or more sites differ in energy after the switching process, a correlated 
electron motion is initiated in the chain. For simplicity, we detune only the energy of the first 
site which leads to the perturbation, 

#i =eo6(t) ™o,a • (6) 

<T=U 

For t > and arbitrary eo, the perturbed system Hq + Hi will initially show a depopulation 
of the first site followed by an accumulation of density on the second. Subsequently, also the 
density on the remaining sites will change with time, and, finally, all occupations will start to 
oscillate. In the case of eo -C 1, the population change will be small, such that the dynamics 
should be well characterized by the linear response properties of the chain. For e > 1, however, 
we expect nonlinear effects to be of importance. In the following, we consider the case eo = 5.0 
and resort to the zero-temperature limit f3 — > oo. Initially, at t = 0, the system is prepared in 
the ground state of Hq. 

We measure the response of the Hubbard chain with respect to the time-dependent electron 
density on the first site (s = 0), 

7 CT (t) = (n^)(t)-i f di(nUt)), (7) 
1 Jo 

-f a (io)=[ dt exp(-iwt) 7 <T (t) , (8) 
Jo 



9. 
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(t, t') = TMc{±[t - If}) exp -i / dt h^ F (t) 



where (nfyft) = — ig^Q (t, t). Due to the spin symmetry obeyed by Eq. ([6j), it is 7(0;) = 7^(w) = 
7^(cj). Furthermore, the time t* indicates a finite propagation time used in the numerics and 
is chosen sufficiently long such that it only affects the basic width of the peaks in j(oj) but not 
their position. 



3.1. Exact dynamics 

In order to benchmark the approximate solution of the KBEs under the use of the GKBA and the 
second Born approximation, we compare to exact reference data computable for small L. The 
exact electron dynamics are obtained from propagating the one-particle nonequilibrium Green's 
function independently of the KBEs. To this end, we rewrite the average in the definition ([!]) 
of the NEGF as, 

L 

9sAt,t') = -^- E^^^^^^t^^^l^cc^W^^Ol^t,^^)^ (9) 

7V t ,A^=0 n 



where Z = N Yl n e ^ Elf t' N ±'' n ^ is the partition function, N = + N±, and n labels 
the states in the iV-particle subspace. For small system sizes, the action of the time-dependent 
creation and annihilation operators, 

c^l(t) = U(0,t)c^U(t,0) , (10) 

with, 

rt' 



U(t',t) = exp ( -ij dt (H (t) - pN) j , (11) 



on the energy eigenstates \N^, n) (with energy E^ t N^ ;n ) can be computed numerically 
exactly. For each particle sector involved in Eq. Q, we work in the trivial site basis which 
involves the states |0), |t), |4) and |t|) and apply the Krylov method [25J together with a 
high-order commutator- free exponential time propagation, see [26j . 

In Figure [T] a) and b), we show exact results based on Eqs. (10) and (11), where the 



consideration of the thermodynamic ground state limits the trace in Eq. (|9j) to iV-j- = N± 
The thick curves show the response for a two-site chain (L = 2) at different repulsive 

interaction strengths U according to Eq. For U = 0, see the black curve in the upper panel, 
there exist two peaks — one at ujq = 5.385 and one at 2ujq = 10.770. Interestingly, for U > 0, 
the energetically lowest peak at ujq develops into a double-peak structure of which the left part 
steadily has a larger spectral weight. For larger U there is even a weight difference of several 
orders of magnitude, cf. Figure [l]b). On the contrary, the peak at 2wo changes only slightly 
with U and almost maintains its spectral weight. 

To understand the spectrum, we recall the switch between two time-independent 
Hamiltonians given by Eqs. ^ and |6]) without a ramp function. In this context, the initial 
state (given by the ground state of ^0) can be expressed as a superposition of eigenstates of 
H2 = Ho + Hi for t > 0. Consequently, the transition frequencies between the eigenstates of 
H2 determine the system's nonequilibrium dynamics, and we expect them therefore to appear 
in j(oj). The inset in Figure [T] a) shows the eigenenergies in the asymmetric chain with eo = 5.0 
and ei = as a function of U, and the thin vertical lines in Figure [T] a) and b) indicate the 
associated excitation frequencies for the selected interaction strengths. Some of the vertical lines 
indeed coincide with the peaks in 7(a)). More precisely, we find three out of six transitions: WAd 
wcd and wad, see the arrows in the inset. In the limit U — > 0, cjac and wcd become degenerate 
which explains the double-peak structure for small, finite U. 
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Figure 1. Exact nonlinear density-response spectra 7(w) = 7^(k>) = 7~K<^>) for a two-site chain 
at eo = 5.0 and half-filling, a) The noninteracting case U = and low to moderate interaction 
strengths U = 0.5, 1.0 and 1.5. b) Moderate towards strong interaction U = 2.0, 3.0, 4.0 and 
5.0. Inset in a): Static energy spectrum of the asymmetric two-site chain with eo = 5.0 and 
e\ = 0. From left to right, the arrows indicate the transition frequencies load, w cd an d wac- 
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Figure 2. a) and b): Nonlinear density-response spectra 7(0;) in Hartree-Fock (HF) 
approximation. The parameters and the thin vertical lines indicating the exact transition 
frequencies are the same as in Figure [TJ The inset in a) compares the [/-dependence of the 
exact transition frequencies u;ac an d ^cd (solid lines) to the Hartree-Fock result (crosses). The 
dashed line is a linear fit, cf. also Figure[5]and Eq. (14). For U = 1.0, the inset in b) shows the 
emergence of "higher harmonics" of the main peak in the HF solution. 



3.2. Approximate dynamics 

In this section, we investigate the density-response 7(0;) in Hartree-Fock as well as in second 
Born approximation propagating the KBEs Q. In the latter case, we in addition apply the 
GKBA of Eq. @. 

The results of Hartree-Fock (HF) calculations for the two-site chain are summarized in 
Figure [2] a) and b) . What immediately becomes evident is the fact that there is no double- 
peak structure developing for finite U > as in Figure [TJ Instead, we observe only a single peak 
below 00 = 6 which moves towards smaller energies as U is increased. The inset in the upper panel 
indicates that the corresponding peak position changes linearly with the on-site interaction U. 
Moreover, its proximity to the lower black solid line referring to oo^c as a function of U suggests 
that the HF approximation covers only the transition A f-> C, compare with Figure [I] a) (inset). 
On top of that, the HF approximation seems not to describe the structure of 7(0;) for frequencies 
around 00 = 11. Here, the observed peaks are obviously "higher harmonics" of the main peak, 
cf. the inset in Figure [2]b). 

Next, we perform GKBA calculations including correlations on the second Born (2B) level. 
Figure [3] shows the result (dots) for U = 0.5 and compares to the HF solution (dashed line) and 
the exact data (solid line). Interestingly, the 2B+GKBA calculation drastically improves the 
HF result in the sense that, instead of a single peak at ljhf = 5.234, there appears a double-peak 
structure with practically vanishing spectral weight around oouf and almost the correct energy 
difference. However, the individual peaks are broadened [28], and in addition there emerge side 
peaks at distances comparable to the double-peak energy spacing Aoo. Furthermore, we obtain 
a similar structure for 7(0;) between 00 = 9 and 12 but with one central peak. Importantly, the 
position of this peak is not a multiple of the frequency of one of the peaks belonging to the 
double-peak at lower frequencies. Hence, it is not a higher harmonics artifact as observed in the 
Hartree-Fock solution (see also Figure [5]) . 

In Figure [4j we show the nonlinear density-response 7(0;) in 2B+GKBA for different 
interaction strengths and investigate the noninteracting limit U — > 0. For practically all values of 
U shown, we are able to identify the double-peak structure which is absent in Hartree-Fock. This 
is the case, though the broadening of the peaks strongly increases with U. In the opposite case, 
for U approaching zero, the double-peak structure correctly merges into a single peak indicating 
the degenerate transition frequencies wac an d ooqd- In the course of this, also the additional 
side peaks (which throughout have frequency offsets of multiples of Aw) vanish. In general, the 
resolution of the side peaks depends on the window function used in the Fourier transform of 
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Figure 3. Dynamics of the two-site chain as in Figures [T] and [2] for the case of U = 0.5. 
Comparison of the exact density response ~f(oo) to the result in Hartree-Fock (HF) and second 
Born (2B+GKBA) approximation. The thin black vertical lines are the same as the red dashed 
lines in Figures [I] a) or [2] a) . 
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Figure 4. Nonlinear density-response spectrum j(u)) for a two-site chain at eo = 5.0 and 
half-filling for selected interaction strengths U in second Born approximation and using the 
GKBA. For the case of U = 1.0, the crosses (and the green line) indicate the result for full, self- 
consistent two-time second Born calculation without applying the GKBA. In the time domain, 
the corresponding function j(t) is known to be strongly damped, compare with [24] . 

Eq. ^ — while in Figure [3] we use a Hamming window, the data shown in Figure [4] are obtained 
by applying a Hann window. However, we note that the characteristic form and broadening 
of the peaks in 2B+GKBA is independent of the window function and also independent of t* , 
cf. Eq. @. 

Finally, in Figure|4j we compare the second Born result for U = 1.0 also to full, self-consistent 
two-time NEGF calculations without the GKBA for /3 = 100 (see the green crosses). Here, 
7(0;) is extremely broad such that we do not obtain precise information about the dynamical 
properties and barely can extract frequencies for the involved transitions. The reason for this 
extreme broadening is a strong damping of the time-dependent occupations (nf )(t) in 2B, for 
details see [24] . 

4. Discussion 

The comparison of Figure [T] with Figures [2] to |4j allows for the statement that 2B+GKBA seems 
to capture the correct physics (in contrast to the Hartree-Fock results and full 2B calculations) 
although there appear additional artifacts in 7(w) such as broadening and side peaks. The main 
result is summarized in the left panel of Figure [5} It shows the exact and approximate transition 
frequencies waCi ^cd and wad as function of the interaction strength U . 

In order to explain the failure of Hartree-Fock in more detail, we investigate the character 
of the excitation processes involved, recall the energy spectrum in the inset of Figure [I] a) and 
see the right panel in Figure |5j For the asymmetric chain with eo = 5.0 and ei = described 
by Hamiltonian H% = Hq + Hi for t > 0, it is easily shown [27J that the excited state D is a 
doubly-excited state relative to the ground state A, whereas the states B and C are singly-excited 
states. With the same argument, state D is furthermore also a doubly-excited state relative to 
the ground state (gs) of Hq. This has the following consequences: 

(i) The energetically largest transition A o D will not appear in the HF solution as double 
excitations are generally not included in any time-dependent Hartree-Fock calculation |21j . 

(ii) The transitions A o C and C f-> D are of one-electron character and therefore in principle 




Figure 5. Left : Comparison of the exact transition frequencies waCj ^cd arid cjad (black solid 
lines) to the results obtained in HF (blue) and 2B+GKBA (red). The thin blue and red lines 
are linear fits to the data points: ojxy(U) = aU + wxy,[/=o> cf. Eqs. (14) and (13). Right: 
Schematic view of the excitation dynamics starting from the ground state (gs) of Hamiltonian 
Hq. For t > 0, we switch to the new Hamiltonian H2 = Hq + H\ , cf. Eqs. \f>\ and (6k. 



describable by Hartree-Fock. The reason why, nonetheless, one observes only the transition 
A f> C is that the function 7(0;) combined with the excitation ^ only probes energy 
differences between states that are populated already at time t = 0. According to (i), 
however, the state D is never populated in HF. 

Points (i) and (ii), clearly explain the simplicity and shortcomings of the HF solution for 
the nonlinear response spectra 7(0;). Moreover, an analysis of the exact eigenstates of the two 
Hamiltonians Hq and H2 reveals that there is finite overlap only between the ground state of Hq 
and the states A, C and D of H2 (all being singlets). The vanishing overlap with the eigenstate 
B (the triplet with S z = [29]) is responsible for the fact that we do not observe transitions 
involving state B in the density-response. In terms of the wave function, the dynamics of the 
system for t > is therefore, 

=CA e- i ^|A)+ Cc e- ii?c '|C)+ CD e- i£D '|D) , (12) 

with |X) being the eigenstates of H2 (having energy Ex) and cx = (X|gs) denoting the expansion 
coefficients with respect to the ground state |gs) of Hq. 

In contrast to the HF approximation, a second-order self-energy (as 2B) has the ability 
to describe double excitations, see [22]. For this reason, the 2B+GKBA calculations capture 
the transition cjad and, accordingly, also wcd- In this context, we again emphasize that the 
obtained frequency for the transition A -H- D is not a multiple of any other frequency observed 
in 2B+GKBA, compare the red dots and the blue squares in the left panel of Figure [5] The 
linear fits to the data points are, 

^ = -0.92817 + 5.385, Wqd = 0-483 U + 5.385 , oj 2 ^ = -0.385 U + 10.770 . (13) 



In Hartree-Fock, we obtain, 



o;HF _ _ .355 u + 5.385 , wfg = -0.710 Z7 + 10.770 = 2 wfg . (14) 



With the result of Eq. (13), the second Born approximation plus GKBA gives a complete 
picture of the nonequilibrium dynamics in the considered Hubbard chain. However, there is one 
additional effect in 2B+GKBA we cannot ignore. This concerns the approximate value for the 
double excitation frequency wad which behaves differently as function of U than the values for 
the single excitations wac & n d ojqd '■ While the approximate frequencies wac an d wcd follow the 
exact results as U is increased, we observe a negative slope (—0.385 in the linear fit) as function 
of U for wad- This is opposite to the exact solution which reveals a clearly positive slope and 
minor absolute value of dwAD/df/ for small U < 2, see the upper solid black line in Figure [5] 
(left panel). 



5. Summary 

In conclusion, we have investigated a small Hubbard-type system with respect to its excitation 
properties beyond the linear response regime by propagating the Kadanoff-Baym equations in HF 
and 2B+GKBA. In the course of this, we found clear signals of electron-electron correlations, 
and the absence of damped solutions for the density-response 7(i) in 2B+GKBA allowed us 
to compute the approximate transition frequencies as function of the interaction strength. 
Comparisons with the exact Green's function furthermore showed that 2B+GKBA produces 
reasonable results with only minor deficiencies regarding broadening and side peak effects. 
Finally, our analysis underlined the importance of double excitations in finite systems. 
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